#ifndef __OCPP_Cesium_Ellipsoid_H__
#define __OCPP_Cesium_Ellipsoid_H__

#include "CesiumMath.h"
#include "Cartesian2.h"
#include "Cartesian3.h"
#include "Cartographic.h"
#include "scaleToGeodeticSurface.h"

namespace OCPP
{
	namespace Cesium
	{
		inline void initialize(Ellipsoid& ellipsoid, double x = 0.0, double y = 0.0, double z = 0.0);

		class _CesiumExport Ellipsoid
		{
		public:
			double x = 0.0, y = 0.0, z = 0.0;
			Cartesian3 _radii = Cartesian3();
			Cartesian3 _radiiSquared = Cartesian3();
			Cartesian3 _radiiToTheFourth = Cartesian3();
			Cartesian3 _oneOverRadii = Cartesian3();
			Cartesian3 _oneOverRadiiSquared = Cartesian3();
			double _minimumRadius = 0.0;
			double _maximumRadius = 0.0;
			double _centerToleranceSquared = 0.0;
			double _squaredXOverSquaredZ = 0.0;

			Ellipsoid(double _x = 0.0, double _y = 0.0, double _z = 0.0)
			{
				initialize(*this, _x, _y, _z);
			}

			Ellipsoid(const Ellipsoid& other)
			{
				this->x = other.x;
				this->y = other.y;
				this->z = other.z;
				this->_radii = other._radii;
				this->_radiiSquared = other._radiiSquared;
				this->_radiiToTheFourth = other._radiiToTheFourth;
				this->_oneOverRadii = other._oneOverRadii;
				this->_oneOverRadiiSquared = other._oneOverRadiiSquared;
				this->_minimumRadius = other._minimumRadius;
				this->_maximumRadius = other._maximumRadius;
				this->_centerToleranceSquared = other._centerToleranceSquared;
				this->_squaredXOverSquaredZ = other._squaredXOverSquaredZ;
			}

			Ellipsoid& operator = (const Ellipsoid& other)
			{
				this->x = other.x;
				this->y = other.y;
				this->z = other.z;
				this->_radii = other._radii;
				this->_radiiSquared = other._radiiSquared;
				this->_radiiToTheFourth = other._radiiToTheFourth;
				this->_oneOverRadii = other._oneOverRadii;
				this->_oneOverRadiiSquared = other._oneOverRadiiSquared;
				this->_minimumRadius = other._minimumRadius;
				this->_maximumRadius = other._maximumRadius;
				this->_centerToleranceSquared = other._centerToleranceSquared;
				this->_squaredXOverSquaredZ = other._squaredXOverSquaredZ;

				return (Ellipsoid&)*this;
			}


			/**
			 * Computes an Ellipsoid from a Cartesian specifying the radii in x, y, and z directions.
			 *
			 * @param {Cartesian3} [cartesian=Cartesian3::ZERO] The ellipsoid's radius in the x, y, and z directions.
			 * @param {Ellipsoid} [result] The object onto which to store the result, or undefined if a new
			 *                    instance should be created.
			 * @returns {Ellipsoid} A Ellipsoid instance.
			 *
			 * @exception {DeveloperError} All radii components must be greater than or equal to zero.
			 *
			 * @see Ellipsoid.WGS84
			 * @see Ellipsoid.UNIT_SPHERE
			 */
			static Ellipsoid fromCartesian3(Cartesian3 cartesian) {
				Ellipsoid result;
				initialize(result, cartesian.x, cartesian.y, cartesian.z);
				return result;
			};

			/**
			 * An Ellipsoid instance initialized to the WGS84 standard.
			 *
			 * @type {Ellipsoid}
			 * @constant
			 */
			const static Ellipsoid WGS84;

			/**
			 * An Ellipsoid instance initialized to radii of (1.0, 1.0, 1.0).
			 *
			 * @type {Ellipsoid}
			 * @constant
			 */
			const static Ellipsoid UNIT_SPHERE;

			/**
			 * An Ellipsoid instance initialized to a sphere with the lunar radius.
			 *
			 * @type {Ellipsoid}
			 * @constant
			 */
			const static Ellipsoid MOON;

			/**
			 * The number of elements used to pack the object into an array1.
			 * @type {number}
			 */
			const static int packedLength = Cartesian3::packedLength;

			/**
			 * Stores the provided instance into the provided array1.
			 *
			 * @param {Ellipsoid} value The value to pack.
			 * @param {number[]} array1 The array1 to pack into.
			 * @param {number} [startingIndex=0] The index into the array1 at which to start packing the elements.
			 *
			 * @returns {number[]} The array1 that was packed into
			 */
			static DoubleVector pack(Ellipsoid value, DoubleVector array1, size_t startingIndex = 0) {
				Cartesian3::pack(value._radii, array1, startingIndex);
				return array1;
			};

			/**
			 * Retrieves an instance from a packed array1.
			 *
			 * @param {number[]} array1 The packed array1.
			 * @param {number} [startingIndex=0] The starting index of the element to be unpacked.
			 * @param {Ellipsoid} [result] The object into which to store the result.
			 * @returns {Ellipsoid} The modified result parameter or a Ellipsoid instance if one was not provided.
			 */
			static Ellipsoid unpack(DoubleVector array1, size_t startingIndex = 0) {

				auto radii = Cartesian3::unpack(array1, startingIndex);
				return Ellipsoid::fromCartesian3(radii);
			};

			/**
			 * Computes the unit vector directed from the center of this ellipsoid toward the provided Cartesian position.
			 * @function
			 *
			 * @param {Cartesian3} cartesian The Cartesian for which to to determine the geocentric normal.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided.
			 */
			static Cartesian3 geocentricSurfaceNormal(Cartesian3 cartesian)
			{
				return Cartesian3::normalize(cartesian);
			}

			/**
			 * Computes the normal of the plane tangent to the surface of the ellipsoid at the provided position.
			 *
			 * @param {Cartographic} cartographic The cartographic position for which to to determine the geodetic normal.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided.
			 */
			static Cartesian3 geodeticSurfaceNormalCartographic(
				Cartographic cartographic) {

				auto longitude = cartographic.longitude;
				auto latitude = cartographic.latitude;
				auto cosLatitude = cos(latitude);

				auto x = cosLatitude * cos(longitude);
				auto y = cosLatitude * sin(longitude);
				auto z = sin(latitude);

				Cartesian3 result;
				result.x = x;
				result.y = y;
				result.z = z;
				return Cartesian3::normalize(result);
			};

			/**
			 * Computes the normal of the plane tangent to the surface of the ellipsoid at the provided position.
			 *
			 * @param {Cartesian3} cartesian The Cartesian position for which to to determine the surface normal.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided, or undefined if a normal cannot be found.
			 */
			Cartesian3 geodeticSurfaceNormal(Cartesian3 cartesian) {
				Cartesian3 result;

				if (Cartesian3::equalsEpsilon(cartesian, Cartesian3::ZERO, CesiumMath::EPSILON14)) {
					return result;
				}

				result = Cartesian3::multiplyComponents(
					cartesian,
					_oneOverRadiiSquared
				);
				return Cartesian3::normalize(result);
			};

			/**
			 * Converts the provided cartographic to Cartesian representation.
			 *
			 * @param {Cartographic} cartographic The cartographic position.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided.
			 *
			 * @example
			 * //Create a Cartographic and determine it's Cartesian representation on a WGS84 ellipsoid.
			 * const position = new Cesium.Cartographic(Cesium.Math.toRadians(21), Cesium.Math.toRadians(78), 5000);
			 * const cartesianPosition = Cesium.Ellipsoid.WGS84.cartographicToCartesian(position);
			 */
			Cartesian3 cartographicToCartesian(Cartographic cartographic) {
				//`cartographic is required` is thrown from geodeticSurfaceNormalCartographic.
				auto n = geodeticSurfaceNormalCartographic(cartographic);
				auto k = Cartesian3::multiplyComponents(_radiiSquared, n);
				auto gamma = sqrt(Cartesian3::dot(n, k));
				k = Cartesian3::divideByScalar(k, gamma);
				n = Cartesian3::multiplyByScalar(n, cartographic.height);

				return Cartesian3::add(k, n);
			};

			/**
			 * Converts the provided array1 of cartographics to an array1 of Cartesians.
			 *
			 * @param {Cartographic[]} cartographics An array1 of cartographic positions.
			 * @param {Cartesian3[]} [result] The object onto which to store the result.
			 * @returns {Cartesian3[]} The modified result parameter or a new Array instance if none was provided.
			 *
			 * @example
			 * //Convert an array1 of Cartographics and determine their Cartesian representation on a WGS84 ellipsoid.
			 * const positions = [new Cesium.Cartographic(Cesium.Math.toRadians(21), Cesium.Math.toRadians(78), 0),
			 *                  new Cesium.Cartographic(Cesium.Math.toRadians(21.321), Cesium.Math.toRadians(78.123), 100),
			 *                  new Cesium.Cartographic(Cesium.Math.toRadians(21.645), Cesium.Math.toRadians(78.456), 250)];
			 * const cartesianPositions = Cesium.Ellipsoid.WGS84.cartographicArrayToCartesianArray(positions);
			 */
			std::vector<Cartesian3> cartographicArrayToCartesianArray(
				std::vector<Cartographic> cartographics
			) {
				std::vector<Cartesian3> result;
				auto length = cartographics.size();
				result.resize(length);

				for (auto i = 0; i < length; i++) {
					result[i] = cartographicToCartesian(cartographics[i]);
				}
				return result;
			};

			/**
			 * Converts the provided cartesian to cartographic representation.
			 * The cartesian is undefined at the center of the ellipsoid.
			 *
			 * @param {Cartesian3} cartesian The Cartesian position to convert to cartographic representation.
			 * @param {Cartographic} [result] The object onto which to store the result.
			 * @returns {Cartographic} The modified result parameter, new Cartographic instance if none was provided, or undefined if the cartesian is at the center of the ellipsoid.
			 *
			 * @example
			 * //Create a Cartesian and determine it's Cartographic representation on a WGS84 ellipsoid.
			 * const position = new Cesium.Cartesian3(17832.12, 83234.52, 952313.73);
			 * const cartographicPosition = Cesium.Ellipsoid.WGS84.cartesianToCartographic(position);
			 */
			Cartographic cartesianToCartographic(Cartesian3 cartesian) {
				//`cartesian is required.` is thrown from scaleToGeodeticSurface
				auto p = scaleToGeodeticSurface(cartesian);

				auto n = geodeticSurfaceNormal(p);
				auto h = Cartesian3::subtract(cartesian, p);

				auto longitude = atan2(n.y, n.x);
				auto latitude = asin(n.z);
				auto height =
					CesiumMath::sign(Cartesian3::dot(h, cartesian)) * Cartesian3::magnitude(h);

				Cartographic result;
				result.longitude = longitude;
				result.latitude = latitude;
				result.height = height;
				return result;
			};

			/**
			 * Converts the provided array1 of cartesians to an array1 of cartographics.
			 *
			 * @param {Cartesian3[]} cartesians An array1 of Cartesian positions.
			 * @param {Cartographic[]} [result] The object onto which to store the result.
			 * @returns {Cartographic[]} The modified result parameter or a new Array instance if none was provided.
			 *
			 * @example
			 * //Create an array1 of Cartesians and determine their Cartographic representation on a WGS84 ellipsoid.
			 * const positions = [new Cesium.Cartesian3(17832.12, 83234.52, 952313.73),
			 *                  new Cesium.Cartesian3(17832.13, 83234.53, 952313.73),
			 *                  new Cesium.Cartesian3(17832.14, 83234.54, 952313.73)]
			 * const cartographicPositions = Cesium.Ellipsoid.WGS84.cartesianArrayToCartographicArray(positions);
			 */
			std::vector<Cartographic> cartesianArrayToCartographicArray(
				std::vector<Cartesian3> cartesians
			) {
				auto length = cartesians.size();
				std::vector<Cartographic> result;
				result.resize(length);
				for (auto i = 0; i < length; ++i) {
					result[i] = cartesianToCartographic(cartesians[i]);
				}
				return result;
			};

			/**
			 * Scales the provided Cartesian position along the geodetic surface normal
			 * so that it is on the surface of this ellipsoid.  If the position is
			 * at the center of the ellipsoid, this function returns undefined.
			 *
			 * @param {Cartesian3} cartesian The Cartesian position to scale.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter, a new Cartesian3 instance if none was provided, or undefined if the position is at the center.
			 */
			Cartesian3 scaleToGeodeticSurface(Cartesian3 cartesian) {
				return scaleToGeodeticSurface1(
					cartesian,
					_oneOverRadii,
					_oneOverRadiiSquared,
					_centerToleranceSquared
				);
			};

			/**
			 * Scales the provided Cartesian position along the geocentric surface normal
			 * so that it is on the surface of this ellipsoid.
			 *
			 * @param {Cartesian3} cartesian The Cartesian position to scale.
			 * @param {Cartesian3} [result] The object onto which to store the result.
			 * @returns {Cartesian3} The modified result parameter or a new Cartesian3 instance if none was provided.
			 */
			Cartesian3 scaleToGeocentricSurface(Cartesian3 cartesian) {
				Cartesian3 result;

				auto positionX = cartesian.x;
				auto positionY = cartesian.y;
				auto positionZ = cartesian.z;
				Cartesian3 oneOverRadiiSquared = _oneOverRadiiSquared;

				auto beta =
					1.0 /
					sqrt(
						positionX * positionX * oneOverRadiiSquared.x +
						positionY * positionY * oneOverRadiiSquared.y +
						positionZ * positionZ * oneOverRadiiSquared.z
					);

				return Cartesian3::multiplyByScalar(cartesian, beta);
			};

			/**
			 * Transforms a Cartesian X, Y, Z position to the ellipsoid-scaled space by multiplying
			 * its components by the result of {@link Ellipsoid#oneOverRadii}.
			 *
			 * @param {Cartesian3} position The position to transform.
			 * @param {Cartesian3} [result] The position to which to copy the result, or undefined to create and
			 *        return a new instance.
			 * @returns {Cartesian3} The position expressed in the scaled space.  The returned instance is the
			 *          one passed as the result parameter if it is not undefined, or a new instance of it is.
			 */
			Cartesian3 transformPositionToScaledSpace(
				Cartesian3 position
			) {
				Cartesian3 result;
				return Cartesian3::multiplyComponents(position, _oneOverRadii);
			};

			/**
			 * Transforms a Cartesian X, Y, Z position from the ellipsoid-scaled space by multiplying
			 * its components by the result of {@link Ellipsoid#radii}.
			 *
			 * @param {Cartesian3} position The position to transform.
			 * @param {Cartesian3} [result] The position to which to copy the result, or undefined to create and
			 *        return a new instance.
			 * @returns {Cartesian3} The position expressed in the unscaled space.  The returned instance is the
			 *          one passed as the result parameter if it is not undefined, or a new instance of it is.
			 */
			Cartesian3 transformPositionFromScaledSpace(
				Cartesian3 position
			) {
				Cartesian3 result;
				return Cartesian3::multiplyComponents(position, _radii);
			};

			/**
			 * Compares this Ellipsoid against the provided Ellipsoid componentwise and returns
			 * <code>true</code> if they are equal, <code>false</code> otherwise.
			 *
			 * @param {Ellipsoid} [right] The other Ellipsoid.
			 * @returns {boolean} <code>true</code> if they are equal, <code>false</code> otherwise.
			 */
			bool equals(Ellipsoid right) {
				return (Cartesian3::equals(_radii, right._radii)
					);
			};

			/**
			 * Creates a string representing this Ellipsoid in the format '(radii.x, radii.y, radii.z)'.
			 *
			 * @returns {string} A string representing this ellipsoid in the format '(radii.x, radii.y, radii.z)'.
			 */
			String toString() {
				return _radii.toString();
			};

			/**
			 * Computes a point which is the intersection of the surface normal with the z-axis.
			 *
			 * @param {Cartesian3} position the position. must be on the surface of the ellipsoid.
			 * @param {number} [buffer = 0.0] A buffer to subtract from the ellipsoid size when checking if the point is inside the ellipsoid.
			 *                                In earth case, with common earth datums, there is no need for this buffer since the intersection point is always (relatively) very close to the center.
			 *                                In WGS84 datum, intersection point is at max z = +-42841.31151331382 (0.673% of z-axis).
			 *                                Intersection point could be outside the ellipsoid if the ratio of MajorAxis / AxisOfRotation is bigger than the square root of 2
			 * @param {Cartesian3} [result] The cartesian to which to copy the result, or undefined to create and
			 *        return a new instance.
			 * @returns {Cartesian3 | undefined} the intersection point if it's inside the ellipsoid, undefined otherwise
			 *
			 * @exception {DeveloperError} position is required.
			 * @exception {DeveloperError} Ellipsoid must be an ellipsoid of revolution (radii.x == radii.y).
			 * @exception {DeveloperError} Ellipsoid.radii.z must be greater than 0.
			 */
			Cartesian3 getSurfaceNormalIntersectionWithZAxis(
				Cartesian3 position,
				double buffer = 0.0
			) {
				Cartesian3 result;

				if (
					!CesiumMath::equalsEpsilon(
						_radii.x,
						_radii.y,
						CesiumMath::EPSILON15
					)
					) {
					/*				throw new DeveloperError(
										"Ellipsoid must be an ellipsoid of revolution (radii.x == radii.y)"
									);*/
				}

				auto squaredXOverSquaredZ = _squaredXOverSquaredZ;

				result.x = 0.0;
				result.y = 0.0;
				result.z = position.z * (1 - squaredXOverSquaredZ);

				if (std::abs(result.z) >= _radii.z - buffer) {
					return result;
				}

				return result;
			};

			/**
			 * Computes the ellipsoid curvatures at a given position on the surface.
			 *
			 * @param {Cartesian3} surfacePosition The position on the ellipsoid surface where curvatures will be calculated.
			 * @param {Cartesian2} [result] The cartesian to which to copy the result, or undefined to create and return a new instance.
			 * @returns {Cartesian2} The local curvature of the ellipsoid surface at the provided position, in east and north directions.
			 *
			 * @exception {DeveloperError} position is required.
			 */
			Cartesian2 getLocalCurvature(Cartesian3 surfacePosition);

			inline const static DoubleVector abscissas = {
				0.14887433898163,
					0.43339539412925,
					0.67940956829902,
					0.86506336668898,
					0.97390652851717,
					0.0
			};
			inline const static DoubleVector  weights = {
				0.29552422471475,
				0.26926671930999,
				0.21908636251598,
				0.14945134915058,
				0.066671344308684,
				0.0,
			};

			using BasicCallback = std::function<double(double)>;

			/**
			 * Compute the 10th order Gauss-Legendre Quadrature of the given definite integral.
			 *
			 * @param {number} a The lower bound for the integration.
			 * @param {number} b The upper bound for the integration.
			 * @param {Ellipsoid~RealValuedScalarFunction} func The function to integrate.
			 * @returns {number} The value of the integral of the given function over the given domain.
			 *
			 * @private
			 */
			double gaussLegendreQuadrature(double a, double b, Ellipsoid::BasicCallback func) {
				// The range is half of the normal range since the five weights add to one (ten weights add to two).
				// The values of the abscissas are multiplied by two to account for 
				auto xMean = 0.5 * (b + a);
				auto xRange = 0.5 * (b - a);

				auto sum = 0.0;
				for (auto i = 0; i < 5; i++) {
					auto dx = xRange * abscissas[i];
					sum += weights[i] * (func(xMean + dx) + func(xMean - dx));
				}
				// Scale the sum to the range of x.
				sum *= xRange;
				return sum;
			}

			/**
			 * A real valued scalar function.
			 * @callback Ellipsoid~RealValuedScalarFunction
			 *
			 * @param {number} x The value used to evaluate the function.
			 * @returns {number} The value of the function at x.
			 *
			 * @private
			 */

			 /**
			  * Computes an approximation of the surface area of a rectangle on the surface of an ellipsoid using
			  * Gauss-Legendre 10th order quadrature.
			  *
			  * @param {Rectangle} rectangle The rectangle used for computing the surface area.
			  * @returns {number} The approximate area of the rectangle on the surface of this ellipsoid.
			  */
			  //static double surfaceArea(Rectangle rectangle) {
			  //	auto minLongitude = rectangle.west;
			  //	auto maxLongitude = rectangle.east;
			  //	auto minLatitude = rectangle.south;
			  //	auto maxLatitude = rectangle.north;

			  //	while (maxLongitude < minLongitude) {
			  //		maxLongitude += CesiumMath::TWO_PI;
			  //	}

			  //	auto radiiSquared = _radiiSquared;
			  //	auto a2 = radiiSquared.x;
			  //	auto b2 = radiiSquared.y;
			  //	auto c2 = radiiSquared.z;
			  //	auto a2b2 = a2 * b2;
			  //	return gaussLegendreQuadrature(minLatitude, maxLatitude, function(lat) {
			  //		// phi represents the angle measured from the north pole
			  //		// sin(phi) = sin(pi / 2 - lat) = cos(lat), cos(phi) is similar
			  //		auto sinPhi = cos(lat);
			  //		auto cosPhi = sin(lat);
			  //		return (
			  //			cos(lat) *
			  //			gaussLegendreQuadrature(minLongitude, maxLongitude, function(lon) {
			  //			auto cosTheta = cos(lon);
			  //			auto sinTheta = sin(lon);
			  //			return sqrt(
			  //				a2b2 * cosPhi * cosPhi +
			  //				c2 *
			  //				(b2 * cosTheta * cosTheta + a2 * sinTheta * sinTheta) *
			  //				sinPhi *
			  //				sinPhi
			  //			);
			  //		})
			  //			);
			  //	});
			  //};

			  /**
			   * Compares the provided Cartesians componentwise and returns
			   * <code>true</code> if they are equal, <code>false</code> otherwise.
			   *
			   * @param {Ellipsoid} [left] The first Cartesian.
			   * @param {Ellipsoid} [right] The second Cartesian.
			   * @returns {boolean} <code>true</code> if left and right are equal, <code>false</code> otherwise.
			   */
			static bool equals(Ellipsoid& left, Ellipsoid& right) {
				Ellipsoid result;
				return (
					left.x == right.x &&
					left.y == right.y
					);
			};

			/**
			 * Compares the provided Cartesians componentwise and returns
			 * <code>true</code> if they pass an absolute or relative tolerance test,
			 * <code>false</code> otherwise.
			 *
			 * @param {Ellipsoid} [left] The first Cartesian.
			 * @param {Ellipsoid} [right] The second Cartesian.
			 * @param {number} [relativeEpsilon=0] The relative epsilon tolerance to use for equality testing.
			 * @param {number} [absoluteEpsilon=relativeEpsilon] The absolute epsilon tolerance to use for equality testing.
			 * @returns {boolean} <code>true</code> if left and right are within the provided epsilon, <code>false</code> otherwise.
			 */
			bool equalsEpsilon(
				Ellipsoid& left,
				Ellipsoid& right,
				double relativeEpsilon = 0.0,
				double absoluteEpsilon = 0.0
			) {
				return (
					CesiumMath::equalsEpsilon(
						left.x,
						right.x,
						relativeEpsilon,
						absoluteEpsilon
					) &&
					CesiumMath::equalsEpsilon(
						left.y,
						right.y,
						relativeEpsilon,
						absoluteEpsilon
					)
					);
			};

			/**
			 * Compares this Cartesian against the provided Cartesian componentwise and returns
			 * <code>true</code> if they pass an absolute or relative tolerance test,
			 * <code>false</code> otherwise.
			 *
			 * @param {Ellipsoid} [right] The right hand side Cartesian.
			 * @param {number} [relativeEpsilon=0] The relative epsilon tolerance to use for equality testing.
			 * @param {number} [absoluteEpsilon=relativeEpsilon] The absolute epsilon tolerance to use for equality testing.
			 * @returns {boolean} <code>true</code> if they are within the provided epsilon, <code>false</code> otherwise.
			 */
			bool equalsEpsilon(
				Ellipsoid& right,
				double relativeEpsilon,
				double absoluteEpsilon
			) {
				return Ellipsoid::equalsEpsilon(
					*this,
					right,
					relativeEpsilon,
					absoluteEpsilon
				);
			};
		};

		inline void initialize(Ellipsoid& ellipsoid, double x, double y, double z) {
			ellipsoid._radii = Cartesian3(x, y, z);

			ellipsoid._radiiSquared = Cartesian3(x * x, y * y, z * z);

			ellipsoid._radiiToTheFourth = Cartesian3(
				x * x * x * x,
				y * y * y * y,
				z * z * z * z
			);

			ellipsoid._oneOverRadii = Cartesian3(
				x == 0.0 ? 0.0 : 1.0 / x,
				y == 0.0 ? 0.0 : 1.0 / y,
				z == 0.0 ? 0.0 : 1.0 / z
			);

			ellipsoid._oneOverRadiiSquared = Cartesian3(
				x == 0.0 ? 0.0 : 1.0 / (x * x),
				y == 0.0 ? 0.0 : 1.0 / (y * y),
				z == 0.0 ? 0.0 : 1.0 / (z * z)
			);

			ellipsoid._minimumRadius = (std::min)(x, (std::min)(y, z));

			ellipsoid._maximumRadius = (std::max)(x, (std::max)(y, z));

			ellipsoid._centerToleranceSquared = CesiumMath::EPSILON1;

			if (ellipsoid._radiiSquared.z != 0) {
				ellipsoid._squaredXOverSquaredZ =
					ellipsoid._radiiSquared.x / ellipsoid._radiiSquared.z;
			}
		}

	}
}

#endif

